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ABSTRACT 

A  system  of  non-linear  partial  differential  equations 
for  the  neutron  density  and  temperature  in  a  circulating 
fuel  reactor  is  studied,   A  steady-state  solution  of  these 
equations  is  obtained  together  with  conditions  which  guaran- 
tee its  existence  and  uniqueness.   Further  conditions  are 
obtained  which  are  sufficient  to  insure  the  boundedness  of 
the  temperature  for  all  time.  Explicit  bounds  are  obtained 
when  the  excess  multiplication  factor  is  a  linear  function 
of  temperature, 

A  numerical  method  for  the  solution  of  the  non-linear 
equations  is  presented.   This  procedure  is  shown  to  be 
stable  and  convergent  under  stated  conditions  on  the  mesh» 
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Section  1.    Introduction 

A  system  of  non-linear  partial  differential  equations 
for  the  thermal  neutron  density  and  temperature  in  a  circu- 
lating fuel  reactor  has  been  derived  and  studied  by 
Dr.  J,  F]e  ck  [1]  of  the  Brookhaven  National  Laboratory.   In 
the  derivation,  which  is  suminarized  in  section  2,  the  excess 
multiplication  factor,  k  ,  is  taken  to  be  a  linear  function 
of  the  temperature  V,   The  present  report  considers  a  gener- 
al, unspecified,  dependence  of  k   on  temperature. 

In  section  3  sufficient  conditions  are  obtained  on  the 
function  k   (V)  to  guarantee  the  existence  of  non-trivial 
steady  state  solutions.   Additional  conditions  are  stated 
which  are  sufficient  to  insure  that  a  non-constant  steady 
state  is  unique.   The  steady  state  temperature  and  neutron 
density  are  given  implicitly  in  terms  of  integrals  of  k  „(V), 
These  results  are  applied,  as  an  example,  to  the  special  case 
assumed  in  [1].   Here  a  bound  on  the  steady  state  exit  temper- 
ature is  obtained  in  terms  of  the  input  temperature. 

In  section  l\.   the  question  of  boundedness  of  the  tem-per- 
ature  for  all  time  is  considered,   A  relatively  weak  condition 
on  k   (V)  is  shown  to  be  sufficient  to  insure  that  the  temper- 
ature  remain  between  specified  bounds  which  depend  on  the 
initial  and  boundary  data  and  the  form  of  the  function  k   (V), 

6  A 

In  the  case  of  the  linear  temperature  dependence,  assumed  in 
[1],  which  always  satisfies  the  required  condition  for 
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boundedness,  explicit  boiinds  are  obtained, 

A  method  for  the  numerical  computation  of  the  temper- 
ature and  neutron  density  is  presented  in  Section  5»   The 
conditions  for  stability  and  an  estimate  of  the  error  in 
this  scheme  are  examined  in  Section  6,   A  discussion  of  some 
of  the  numerical  results  which  have  been  obtained  with  the 
use  of  this  method  is  to  be  included  in  [1]. 

Section  2,    Formulation 

Let  a  cylindrical  circulating  fuel  reactor  with  axis 
parallel  to  the  X-axis  occupy  the  range  0  <  X  <  A.   The 
thermal  neutron  density  may  be  taken  as  U(T,  X)R(y,  z), 
where  R  satisfies 

R   +  R   +  (B^-6^)R  =  0   , 
yy     ZZ    '    ^  '        * 

(2.1) 

R(y,  z)  =  0  ,  for   (y,  z)  on  the  cylindrical  surface; 

and  U(T,  X)  satisfies 

iJU^  =  h^(U^  +  P^U)  +  [k^^^-l]U  SOI  ,. 

(2.2) 

U(T,  0)  =  U(T,  A)  =  0   . 

Here  we  have  used  the  notation 

IT 

^   =  thermal  neutron  lifetime       »  P  =  /A  » 

L  =  thermal  diffusion  length       ,  ^  =  '^o/l+L^B        , 

(2.3)k  =  infinite  reactor  multiplication,  h^  =  L^/l-fL  B^      , 

— B  l'    ?  2 

"7^  =  Fermi  age  of  thermal  neutrons   ,  k^^  =  ke"   /l+L  B  , 


2 

B  =  reactor  buckling. 
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If  the  fuel  is  flowing  with  constant  velocity  C  >  o  in  the 
positive  X-direction,  its  temperature,  given  by  V(T,  X)E(y,  z}, 
must  satisfy 

Vj  +  CV^  =  U 

V(T,  o)  =  V^ 

Here  V  is  the  constant  input  temperature  and  the  units  are 
o 

appropriately  chosen  such  that  temperature  is  measured  in  C, 
It  is  assumed  that  the  significant  effects  of  temperature  on 
the  reactor  kinetics  may  be  obtained  by  taking  the  excess 
multiplication  factor,  tk  f^-l]  =  k  ,  to  be  a  given  fvmction 
of  temperatiire;  say 

(2.5)  tkg^^-1]  =  e-6V  , 

with   e,    6   >  o.      Equations    (2,2),     (2,14.)    and    (2.5)    i^ay  be  written 
in  dimensionless   form  as 

a)  u.    =  u       +    [l-v]u 

(2.6)  ^  ^^ 

b)  v^  +   cv^  =  u  , 

with  the  boundary  conditions 

a)  u(t,  o)  =  u(t,  a)  =  0 


(2.7) 


b)  v(t,  o)  =  Vq 


The  dimensionless  variables  are  defined  by 

t  =  P'/^  I      ,        X  =   P4  X      .  c  =    ^/p^  0 

(2,8)  Of.  c-  p  p    p 

u=±£u  ,        v  =  °/2V,  where  p^  =   e  +  p'^h'^      , 

.5- 


and  a  and  v  are  ths  dimen^lonless  length  and  input  temper- 
a':ure  respectively.   Equations  (2,1)  -  (2,8)  are  a  suinmary 
o.^   the  formulation  derived  in  [1]-: 

The  temperature  dependence  assumed  in  '2.5)  simplifies 
analytical  attempts  to  investigate  the  resulting-  system  of 
eqiAations,   It  also  assures  a  negative  temperature  coef- 
ficient of  reactivity  for  all  positive  temiperatures  which,  at 
least  intuitively,  seems  sufficient  to  prevent  unbounded 
temperatures  and  neutron  densities.   That  this  is  so  for  the 
temperature,  is  shown  as  an  example  in  Section  Ij.,   However 
it  is  important  to  know  what  general  types  of  dependence  of 
k   on  V  lead  to  bounded  temperatures.   Thus  we  take 

(2-9)    :k^^f-l]  =  k^^^(V)    , 

where  the  functional  form  is  unspecified  and  we  are  not 
restricted  to  negative  temperature  coefficients.   Introducing 
variables  analogous  to  those  of  (2,8),  but  not  necessarily 
dimensionless,  equations  (2.2),  (2,14.)  and  (2,9)  become 
a)     u^  =  u^^^  +  A(v)u 


(2,10) 

^'     ^t   -X 
with  the  boundary  conditions  (2.7).   The  restrictions  imposed 


b)     v^  -'-  cv^  =  u     , 


on  k   (V)  are  given  in  Sections  3  and  [j.  as  conditions  on  A(v), 
Thoy  are  shown  to  be  satisfied  for  A(v)  =  1  -v. 

To  ccm^, lece  the  formulation  v/e  specify  the  initial  con- 
ditions 

a)  L.:o,  x)  =  f(x) 

(2.11) 

b)  v(o,  x)  ■■■=   g(x)     , 


'I'.-/ 


where,    since  u  is   a  density,   we   require 

(2.12)  f (x)    >  o      ,      o  <  X  <  a 

In  order  to  avoid  the  consideration  of  discontinuities  we 
also  require 

(2.13)  g(o)  =  v^   , 

and  assume  that  all  functions  introduced  have  continuous 
derivatives  of  at  least  fourth  order. 

Section  3«    The  Steady  State 

Let  the  steady  state  be  described  by  the  variables 

u(x)  =  lira    u(t,  x) 
t  — >  00 

(3.1) 

v(x)  =  lira    v(t,  x)    , 
t  — >oo 

This  notation  will  be  used  only  in  the  present  section  and 

to  avoid  confusion  the  argimients  (t,  x)  shall  be  included 

for  all  non-steady-state  quantities.   With  the  assumption 

that  lira    u+.(t,  x)  =  lira    v,  (t,  x)  =  0,  we  obtain 
t  — 9>00  ^         t  — >co  ^ 

frora  (2,10) 

a)  u"  +  A(v)u  =  0 

b)  cv'  =  u   • 

The  boundary  conditions  are,  from  (2,7) 

a)  u(o)  =  u(a)  =  0 

(3.3) 

b)  v(o)  =  v^     , 
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In  Section  [j.  it   is   shown  that  u(t,    x)   >  o  for  all   t  >  o. 
Anticipating  this    result  we   require 

(3,[(.)  u(x)    >o     ,      o<x<a        • 

One  solution  of  the  above  equations  is  the  constant  steady 

state:  u  =  o.  v  =  v  ,  We  consider  this  to  be  a  trivial 
'      o 

solution  and  seek  other  non-trivial  steady  states. 

Some  qualitative  features  of  non-trivial  solutions  are 
immediate  consequences  of  (3,2b),  (3»3a)  and  (3.ij-).   They 
imply  that  v(x)  must  be  a  monotonic  increasing  function  of  x 
with  horizontal  slope  at  x  =  o  and  x  =  a.   For  a  complete 
analysis  we  introduce  the  quantities 
(3«5)    V.  =  v(a),  the  exit  temperature; 


(3.6) 


B 


<^'  ^q)  ^  /  ^^^^^^         ' 


(3.7)    G(v,  v^)  =  /  B(^,  v^)d4   ; 


e  o 


-il 


(3.8)I(v;  v^,  v^)  =  -^ 


"f 


-1 


(^.v^)C(v^,  v^).(v^-v^)C(4,  v^) 


/2 


dC 


Eliminating  u  in  (3.2),  integrating  the  resulting  differential 
equation,  and  applying  the  boundary  conditions  (3.3 )>  lead  to 
the  following  result:  A  non-trivial  solution  of  {'^,2)    -  (3.^) 
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is  given  implicitly  by 


X  =  I(v;  v^,  v^) 


(3.9) 


^-1 


u  =  c 


e  o 


/2 


~l 


(v.v^)C(v^,  v^)-(v^-v^)C(v,  v^) 


V2 


provided  there  exists  a  v  >  v  such  that 

e    o 

(3.10)   a  =  I(v^;  Vg,  v^) 


and 
(3.11) 


C(4,  V  )   C(v  ,  V  ) 

2-  <  H 2_  for  V  <  £  <  V 

o  —  ^  —  e 


4-v, 


V  -V 

e  o 


The  constant  solution,  u  =  0,  v  =  v  ,  is  unique  if  there  does 
not  exist  a  number  v  >  v  such  that  (3.10)  and  (3.11)  are 
satisfied.   If  there  exists  only  one  such  number,  then  there 
is  only  one  non-trivial  solution,  and  it  is  given  by  (3»9). 

A  sufficient  condition  for  the  uniqueness  of  a  non- 
constant  steady  state  is  that  I(v;  v,  v  )  be  an  increasing 

function  of  v  in  the  range  v  >  v_.  For  such  cases  condition 

^       o 

(3,10)  requires  that 


(3.12)   a  >  lim 
v 


^^^e'  ^e'  ^^   • 


e  — ►  o 

This  limit  is  obtained  by  introducing 


(3.13) 


e  =  v-^  -  v_ 
&         o 


?  =  eT\  +  v^ 


into   (3.8)   to  obtain 
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I(v^.'  ^e'    %)  =  ^^^c  ■*"  ^'  ^o  ■"  ^>  ^o^ 


(3.11+)  ^ 


r 


A  Taylor-series  expan-ion  gives 


\C(v^   +  e,  v^)-C(v^  +  e>|,  v^) 


-1 


/2 


2  2  3  3 

C(v  +  eni,  V  )  =  -^— ^  A(v  )  +  --^  A'  (v^  +  Qtri)    ,  o  <  6  <  1  , 

'O      V'O        p        O        /        O        V'»—    —    » 

where  we  have  used  the  definitions  (3«6)  and  (3.7).   Using  this 
result  in  (3ali|)  and  passing  to  the  limit  yields 


Tf  ^  1     r  d^ 

lim     I(v^;  v^,  v^)  =  /   ^ 

e   e   o 


Thus  condition  (3,12),  which  is  necessary  for  the  existence 
of  a  non-trivial  steady  state  in  the  present  case,  becomes 

(3.15)    a  >  V/aT^ 

In  terms  of  the  dimensional  variables  this  condition  is 
simply  kg^(v  )  >  o   .   If  I(v;  v,  v  )  is  a  decreasing  function 
of  V  for  V  >  V  ,  the  above  condition  with  the  inequality  re- 
versed is  necessary  for  uniqueness. 

The  condition  (3. 11)  has  a  simple  geometric  interpre- 
tation.  Consider  the  curve  y  =  C(^,  v  )  for  4  >  v   ,  Then 

v  must  be  a  number  >  v   such  that  this  curve  lies  below  the 
e  o 

line  with  slope  C(v  ,  v  )/(v  -v  ),  passing  through  the  point 

GO       GO 

^  =  V  ,    y  =   0,      Prom  formula    (3.9)    this   slope   can  be   shown 
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to  have  the  Interpretation 

(3.16)   2 2-.=  v"(o)>o  , 

e  o 

The  non-negativeness  results  from  the  fact  that  -^   =  0  at 

4  =  V  ,   Similarly  we  deduce  from  (3,11) 

(3.17)  ^  I       =B(v^,  v^)  >    ^'     '^ 
d4  I  4  =  V       ®   °      V  -V 


e'   o' 
e  e   o 


The  conditions  (3,16)  and  (3,17)  are  helpful  in  attempting 

to  verify  (3,11)  for  a  specific  problem.   The  former  may 

also  be  used  as  a  convergence  criterion  in  numerical  work. 

As  an  important  illustration  of  the  above  results  we 

consider  the  case  A(^)  =  1  -  ,^  ,   Then  from  definitions  (3«6) 

and  (3.7) 

B(v,  v^)  =  I  (v-v^)(2-v^.v) 
(3.18) 

C^^'  ^o^  =  Z   fv-v^)2(3-2v^-v)    , 

Since,  from  (3,16)  and  (3,17),  B(Vg,  v^)  >  o,  and  since  for 
non-trivial  steady  states  v  >  v  ,  we  may  deduce  from  the 
first  equation  of  (3.18) 

(3.19)  v^  <  1 

This  is  a  necessary  condition  for  the  existence  of  a  non- 
constant  steady  state.   Then  condition  (3,11)  is  satisfied 
provided 

(3.20)  v^  <  v^  <  -—     , 
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which  furnishes  a  boxind  on  the  steady-state  exit-temperature » 
The  above  inequality  is  easily  obtained  by  using  (3,18)  and 
(3,19)  in  (3.17).   The  integral  in  equation  (3.10)  becomes 
in  the  present  case 


1  -. 

(3.21)    I(v^;  v^,  v^)  =  y  7|(l-7|)(l-v^-  I  [l-ft\])| 


-Va 


dT, 


where  e  =  v  -v  ,   This  is  obviously  an  increasing  function 
e  o 

of  e  and  hence  of  v  .   Therefore  we  require  from  condition 

e  ^ 

(3.15) 

(3.22)        a  >  V/l^o 

.>^o 


From   (3.21)   we   see   that   I(v^;    v^,    v^)   — >  00    as    v^  —*■- 

6     6     O  0        -^ 

and  hence  a  unique  non-constant  steady  state  exists  for  any 

a  >  ^/vHTv^  and  v^  <  1.   Using  (3. 18)  in  (3.9)  this  solution 
is  obtained  in  terms  of  an  elliptic  function.   For  the  case 
V  =  0,  the  solution  has  been  given  by  J.  Fleck  [Ij, 


Section  k.»        Temperature  Bounds 

Let  u(t,  x),  v(t,  x)  be  solutions  of  (2.10)  satisfying 
the  boundary  conditions  (2.7)  and  initial  conditions  (2,11- 
2.13).   We  first  prove  that  u(t,  x)  >  o  which  justifies  (3»k) 
and  is  used  later.   The  proof  uses  the  following  form  of  the 
weak  maximum  principle  for  parabolic  equations:  Let  w(t,  x) 
be  a  solution  of  the  equation 
(i|,l)    ^t  ^  ^xx  "^  ^^^'  ^^^ 

'<■   This  section  is  an  adaptation,  in  somewhat  less  rigorous  form, 
of  a  forthcoming  paper  [3]« 
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in   the    rectangle 

(k-2.)  Rm:    o<x<a     ,        o<t<T      . 

If  0(t,    x)  <  o  and  w(t,  x)  >  o  in  Rm,  then  w  attains  its 

maximum  on  one  of  the  boundary  segments 

CIf.,3)    X  =  0,  t  =  0,  or  X  =  a 

This  principle  has  long  been  known  in  various  forms;  see 
for  example  M.  Picone  [  ] » 

In  the  rectangle  R,p  of  il4.,2)   we  define 

w(t,  x)  =  e"^''^u(t,  x) 
where 

X  =  X(T)  =  max  A(v(t,  x)),  for  (t,  x)  e  Rrp   . 

Then  w  satisfies  (l+.l)  with  0{t,    x)  =  A(v)  -  X  <  o  in  R^  , 
Prom  the  above  maximum  principle  we  may  conclude  that  w, 
and  therefore  u,  assume  negative  values  in  R^  only  if  u 
assumes  a  negative  minimum  on  one  of  the  segments  ([|.«3)t 
But  u,  and  therefore  w,  are  >  o  on  (U.3)  by  the  conditions 
(2.7a),  (2,11a)  and  (2.12).   Therefore  u  >  o  in  Rn,  and  since 
T  is  arbitrary  u(t,  x)  >  o  for  all  t  >  o  , 

To  obtain  bounds  for  v(t,  x)  we  consider  the  rectangle 
Rrp  as  the  sura  of  four  boundaries  and  the  interior  as  follows, 

N:   o<x<a   ,    t=T 

E;x=a       ,    o<t<T 
(U.>k)  S:   o<x<a   ,    t  =  0 

V/:x=0       ,    o<t<T 

I;   o<x<a   ,    o<t<T   « 
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We  also  define  sets  R„  and  Rm  as  the  intersections  of  R^ 
with  the  half  planes  ct  -  x  <  o  and  ct  -  x  >  o  respectively. 


X  =  ct 


S      a 

Figure  1, 

A  lower  bound  is  found  immediately  by  considering  the 
segments  of  slope   /c  from  each  point  p  of  R^  to  a  corre- 
sponding point  p'  on  V/  or  S,   Since  u(t,  x)  >  o  it  follows 
by  integration  of  (2.10b)  along  these  segments  that  v(p)  >  v(p'  ) 
However,  from  (2,7b)  and  (2.11b)  we  have  v(p«  )  =  v^^  on  W  and 
v(p'  )  >  ci  =  JT^in  g(x)  on  S.   Thus, 


(i|.5) 


v(p)  >  v^  for  peRrp 
v(p)  >  a  for  peRm 


or  since  g(o)  =  v  wa  have 
(4,6)    v(p)  >  a  for  peR^ 

To  obtain  an  upper  bound  on  v  we  eliminate  u  from  (2,10) 
and  get 


ik-^7) 


^   "  jljC-^t  ^  ^xx-^S^^'  ^oO 


=  0 


where  B(v,  v  )  is  defined  by  (3.6),   This  equation  implies  that 


.11^. 


.         /  v^ 


•  1 


I  - 


(k-^Q)  v^  =  "^xx  ■*"  ^^^»  "^o^  "^  P(ct-x) 

for  some  continuous  function  P.   To  obtain  an  upper  bound 
on  V  we  require  an  upper  bound  on  F,   Thus  setting  t  =  0  in 
(ij.,8)  we  get  on  S 


(1^.9) 


P(-x)  =  v^(o,  x)  -  v^^(o,  x)  -  B(v(o,  x),  v^) 

=  f(x)  -  cg»  (x)  -  g"(x)  -  B(g(x),  v^)    , 


where  we  have  used  (2,10b)  and  the  initial  conditions  (2,11). 

If  we  define  y  by 

(il.lO)   r  =  min  [-f(x)  +  cg'(x)  +  g"(x)  +  B(g(x),  v  )] 
o<x<a 

then  (i;.9)  yields 

(I;. 11)   F(ct-x)  <  -  Y  for  (t,  x)eR~ 

For  positive  argioments  we  set  x  =  0  in  (l+.B)  and  obtain  on  W 

F(ct)  =  v^(t,  o)  -  v^^(t,  o)  -  B(v(t,  o),  Vq) 

(1+.12) 

=  -  v^^(t,  o)    . 

Here  we  have  used  (2,7b)  to  get  v. (t,  o)  =  0  and,  with  the 

definition  (3.6),  B(Vq,  v^)  =  0.   Since,  from  (2.7a),  u(t,  o)  =  0 

and  as  above  v. (t,  o)  =  0  we  get  v  (t,  o)  =  0  from  (2.10b). 

However,  (k»5)    requires  that  v  >  v  in  a  neighborhood  of  each 

point  of  V/,  and  thus 

V   (t,  o)  >  o   • 
XX   *   '  — 

This   result   and    (I|..12)    require 

(1+.13)        P(ct-x)   <  o  for    (t,    x)eRj  , 
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which  together  with  (i+.ll)  furnish  the  required  bounds  on 
P  In  Rrp  .. 

Let  P  be  a  point  of  Rm  at  which  v  attains  its  maxi- 
muiti  value.   Then  P  lies  In  one  of  the  continua  E,  S,  W, 
or  W  +  I  .   If  PeW,  then  v(P)  =  v  by  the  boundary  condition 
(2.7b),   This  is,  in  a  sense,  a  trivial  case  as  then 
v(t,  x)  =  v^  foi^  (t,  x)eR^,  by  (k»B) »      If  PeS,  then  by  the 
initial  condition  (2.11b) 

(k»lk)        v(P)  =  P  =  max  g(x)   , 

o<x<a 

Suppose  PeN  +  I;  then  since  v(P)  is  a  maximum  and  the 

set  N  +  I  is  open  with  respect  to  x  and  open  from  below  with 

respect  to  t,  we  must  have 

(I|..l5)    v^(P)  =  0,  v^(P)  <  o,  v^(P)  >  o  . 

Taking  ([4. 0)  at  this  point  P  yields 

B(v(P),  v^)  +  F(P)  >  o   , 

or  from  the  bounds  (I4..II)  and  (I-I-.I3)  on  F 

'y  for  PeRrp 
B(v(P),  V  )  >  )  ,  or  in  general 

o  for  PeRj 


>  min  (y,  o)  for  PeN  +  I  , 


If  there  exists  an  M  >  v  such  that 

—  o 


(i^a?) 


B(M,  v^)  =  min  (y,  o) 

B(v,  V  )  <  min  (y,  o)  for  v  >  M 
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we  could  conclude  from  (14..  16)  that 
(ij..l8)    v(P)  <  M   „ 

Since  B(v^,  v^)  =  0,  and  B(v,  v  )  is  a  continuous  function 
of  V,  (I4..17)  is  implied  by  the  condition 

([)..19)    lim  sup  B(v,  v  )  <  min  (y,  o)  , 
V  — >  00       ° 

We  assume  this  condition  to  hold,  which  is  a  restriction  on 
the  given  function  A(v),  and  hence  ([).»18)  follows. 

Finally,  suppose  PeE,   Since  v(P)  is  the  maximum  value 
of  V,  the  definition  (ii-.ij.)  of  E  leads  to  the  inequalities 

But  v^  +  cv^  =  u  =  0  on  E,  and  c  is  >  o  .   Therefore  the 
above  inequalities  yield 

(ij..20)    v^(P)  =  0,  v^(P)  =  0  . 

From  the  maximum  property  of  v(P)  we  may  now  conclude 

(I4-.2I)    v^(P)  <  o  . 

From  (i;.20)  and  (i|.21)  follow  (lj-.l5)  and,  hence,  the  ine- 
quality v(P)  <  M  ,   Summarizing  the  results  for  all  four 
sets  E,  S,  W,  and  N  +  I  we  find  v(P)  <  max  (p,  M)  regardless 
of  where  the  maximum  value  v(P)  is  attained  in  Rm.  However, 
since  T  is  arbitrary  we  have  obtained  a  bound  for  all  t  >  o  ♦ 
Thus  we  have  proven  the  following: 

Theorem.   Let  u(t,  x),  v(t,  x)  be  a  solution  of  (2,10) 
satisfying  the  boundary  conditions  (2,7)  and  initial  conditions 
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C2.ll   -   2,13)  <•      Let   B(v,    v^)    satisfy  condition    (1^.19).      Then 
u(t,    x)    >  o   and  v(t,    x)    is   bounded   for  all   t  >  o   and  o  <  x,  <  a     ; 
V  satisfies   the   inequality 

(l+,22)        a  <  v(t,    x)   _<  max    (p,    M)      , 
where   a    -  min  g(x),    p   =  max  g(x),    and  P4  is    defined  by    (k.,17). 

As    an  explicit    application  of  the    above   theorem  we   con- 
sider the    case  A(v)    =  1   -   v.      The   corresponding  B(v,    v   ) 

is  given   in    (3«l8)..      V/e   see   that  lira         B(v,    v    )    =  -   oo    so 

V  — 3>  oo 

that    (l4.,19)    is    satisfied  for  all  y   and  v^„      The   nvunber  M  of 
condition   (Ij..>17)    is    easily  found   to  be 

/  Vq  for  Y  >  o,    v^  >  1 

([;.23)        M  =  <    2-v^  ^  for  Y  >  o,    v^  <  1 

(^  1  +   [(1-v^)^-2y]    ^^  for  Y  <  o      , 

Section_5»   tumeric al_Solut  ion 
We  consider  the  equations 

a)  R[u,  v]  =  v.  +  cv  -  u  =  0 

(5.1) 

b)  S[u,  v]  =  u^  -  u^  -  0{n,    v)  =  0 

subject  to  the  boundary  and  initial  conditions 

u(o,    t)    -  u(a,    t)    =  0\  u(x,    o)    =  f(x) 

i$,2)  Vt>o;  (o<x<a 

v(o,    t)    =   Vq  j  v(x,    o)    =  g(x) 

Equations  (5*1)  are  more  general  than  those  of  (2,10),  and  the 

restrictions  (2.12)  and  (2.13)  on  the  initial  data  have  been 

dropped.   We  assiime  the  existence  of  a  unique  solution  to  the 

above  problem.   In  the  present  section  we  associate  a  system 

of  finite-difference  equations  with  (5.1)  and  show  that  they 
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have   a  unique    solution. 

On  the    rectangle   R-,   of    (I|-.2)    we    introduce   the   mesh 

/x^  =  ih 
(S»3)        )  >      where 

and  W  and  J  are  fixed  positive  integers.   Let  U^,  V;  be 
defined  on  the  mesh  as  solutions  of  the  difference  equations 

iS.k) 
a) 

kR,  [uJ,V^]  =  (V^-vj"-^)+lx(vJ"-'--V^'''3-)-kU^"-'-=0           ,  l<i<N 
h  i*  1  --  1   i   '  t-^  i    1-1'    i  ' 

b)  Vi<J<J 

kSj^[uj,V^]  =  (U^-U^"^)-X(uJ^3^-2uJ+U^^^)-^(uJ-^,V^"^)=0,l<i<N-l^ 
where 

(5.5)  ,       fj>i;     o       >o<i<N. 

VJ  =  v^     (  V°  =  g(x^) 

In  (5«U)  we  have  introduced 

(5.6)  [I  =   c^A  ,   X  =  Vh^ 

The  solution  of  these  equations  is  easily  obtained  by  in- 
duction.   Thus  we  assume  U^  ,  Vv   known  for  all  j'  <  j  -  1  » 
From  {S*k-)    a)  we  have 

v^  =  (i-u)v^'"-'-  +  [iV^.'}  -f  kuJ'-'-      ,   1  <  1  <  N 

which  uniquely  determines  the  V^  •   We  write  (5.1|)  b)  as 

XU^^3^-(1+2X)UJ  +  XUJ_^  =  -  U^-^  -  0{UJ"^,  VJ"^),  1  <  i  <  N  -  1  . 

This  is  a  system  of  W  -  1  equations  for  as  many  unknowns,  since 
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U^  =  U,%  =  0  by  iS»B)  .   The  coefficient  matrix  is  non-singular, 
which  Insures  a  unique  solution  of  this  system  of  equations. 
For  computational  purposes  this  matrix  may  be  factored  into  a 
product  of  a  lower  -  by  an  upper  -  triangular  matrix.   The 
solution  is  then  obtained  by  evaluating  a  pair  of  two-term 
recursions. 

If  the  last  term  in  (5»^)  b)  is  taken  at  j,  rather  than 
j  -  1,  the  resulting  implicit  equations  would,  in  general, 
not  be  linear,  Hov/ever,  they  could  be  solved  by  an  iteration 
scheme  the  convergence  of  which  has  been  examined  by  M,E,Rose[i|] 

Section  6,   Convergence  and  Error  Estimate 
Let  u(x,  t),  v(x,  t)  be  the  imique  solution  of  (5»1) 
and  (5.2)  in  Rrp  and  U^,  V^  bo  the  solution  of  (S.k)    and  (5.5) 
on  the  mesh  (5''3)»   ^^e  define  the  error  of  the  numerical 
solution  as 


(6.1) 


EJ  =  U^  -  u(x^,  tj) 
pJ  =  V^  -  v(x.,  t^) 


In  the   present   section  we   shall  prove   that    the    condition 

(6.2)  [i  =  c^/h  <  1 

is  sufficient  to  insure  the  uniform  boundedness  of   |Ev|  and 

|Pv|,   The  bounds  are  such  that  if  h  — >o  and  k  — >o  while 

maintaining  condition  (6,2),  then  the  error  uniformly  approaches 

zero. 

The  above  results  will  follovj  from  the  bomidedness  of 
the  solution,  X^,  Y^   of  the  following  set  of  linear  differ- 
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ence  equations  on  the  mesh  iSrj)^^ 

(6.3) 

a) 

\i^i»    YJ,]  =  r^   ,     1  <  i  <  N 

b)  1  <  J  <  J 

Lj^[xJ,  Y^]    =Jil  ,     1  <  i  <  N-1 

c) 

(XJ,  yJ,  X^,  X°,  Y°)  =  given,  for  1  ,  .  j  <  J  and  o  <  i  <  W. 

Here  R,  is  defined  as  in  iS^h)    a)  and 

(6.1;) 

kLAyJ,  Yh  =  (xJ-x?~^)-\(xJ' ,-2x?+xj'  .)-JrVr^''-b^r-^Yi'-^   . 

hi*   1—1   1   '     1+1   1   1-1   1   1    1   i 

The  condition  (6,2)  is  assumed  to  hold  for  the  mesh  \S'>3)» 

To  prove  the  boundedness  of  the  X.  and  Y^  we  introduce: 
^  11 

(6.5) 

a) 

A  =  maxlaYj,  B  =  max|b^I    ^ 

b) 

C  =  max(|rfl,  |4!,  lX°|,  |yJ|,  jX^i,  |Y^1,  Ix|1  )    ; 

c) 

X,  =  maxdX^I,  ^)   ,   Y.  =  max(JYJ!,  i)        ; 

1  ^  J  ■"■■  i 

d) 

X.  =  max  X.     ,    Y.  =  max  Y.,  « 

Then  solving  (6,3)  a)  for  Y^  ,  taking  absolute  values,  and 
using  (6,5)  b)  and  c)  we  obtain,  since  [a  <  1, 

« 
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By  definition  we  have  Y.  -,  >  E,   and  thus 

Applying  this  inequality  recursively  yields,  with  the  aid 

of  (6.5)  b)  and  d) 

(6.6)    Yj  <  JkXj_^  +  (l+jk)4 

This  furnishes  a  boiond  on  the  Y^  in  terms  of  a  bound  on  the 

Xv   for  j'  <  j  .   To  obtain  the  latter  boxmd  we  write  (6.3)  b)  as 

(6.7) 

(x|-xJ"^)-X(xJ^3^-2X^+X^_^)  =  k(a|"^X^-^+bj"^Y["Ve^),  1  <  i  <  N-1. 

Prom  (6,5)  b)  and  c)  either 

(6.8)    X.  =  ,|  or  X.  =  max    \X^.\    >•  4 
•^  -J    l<i<N-l  ^ 

Assume  the  latter  case  to  hold  and  let  i,  1  <  i  <  N-1  ,  be  a 

subscript  for  which  X.  =  |X^|,   If  X^  >  o  we  have  from  (6,7) 

(x|-x|~^)  <  k(a^"^xJ~^  +  b^'^Y^"^  +  ^l) 

or  using  (6,5)  a)-c), 

(XJI  <  |X|-^|  +  k(A|X^-^|  +  B|Y^-^|  +  ^) 
(6,9) 

<  (l+kA)X._j^  +  kBY,  ^  +  kg 

If  X^  <  o,  (6.7)  yields 

(XJ-X^-^)  >  k(a^-^x|-^  +  b|-M-^  +  xi)      , 
or  on  multiplication  by  -1, 
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Since,  in  this  case,  -Xt  >  o,  we  may  take  absolute  values  In 
the  above  inequality  and  obtain  again  (6,9).   By  definition 
X.  ,  >  ^,  and  thus  (6,9)  yields  in  either  of  the  cases  (6,8) 

X.  <  (l+kA)X.^^  +  kBY,  ^  +  k4  , 
From  (6,6)  in  the  above  inequality  and  since  ^*  ^  >  ^i.p  ^® 
obtain 

(6.10)  X.  <  (l+kA)X._,  -f  kB(j-l)kI'.  ,  +  [l+B+B(j-l)k]k§ 

Since  (6,10)  holds  for  each  j  in  1  <  j  <  J  we  have,  using 

(6.5)  d) 

X.  <  (l+kA+kBt.)X,  ^  +  (1+B+Bt.)k^  ,  t.  =  Jk 
Applying  this  inequality  recursively  and  taking  j  =  J  we 
obtain 

(6.11)  ^j  <  (  [l+k(A+BT)]'^  -f  iitllL^l§IIL-:i[l^.B+BT]  )  £, 

(^  (A+BT)  J 

Prom  (6,5)  d)  and  the  above  result,  (6,6)  implies 

(6.12)  Yj  <  TXj^^  +  (1+T)^   . 

These  are  the  required  uniform  bounds  for  X^  and  Y^  on  the 
mesh  (5.3).   A  useful  form  of  (6,11)  is  obtained  by  setting 
nd  using 
(1+  Vn)^  <  e^  for  n,  a  >  o  ; 


m 

k  =  /j  and  using 


thus 

(A+BT)T     ^(A+BT)T  , 
(6.13)   Ij  <   ^  e        +  ^ =i  (1+B+BT) 

(A-fBT) 
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Prom  the  definition  (6,1)  of  the  error  and  (5»1)  a), 
iS'h)    a)  we  have 

Rj^(e|,  F^)  =  Rj^CU^,  V^)  -R^(u(x^,  tj),  v(x^,  tj)) 

=  R(u(x^,  t.),v(x^,  t.))  -Rj^(u(x^,  t .),  v(x^,t .)  ) 
Prom  the  usual  Taylor  expansions  this  becomes 

f6a/+) 

Rh(E^,Pi)=r.^.=  |v^^(x^,t..e^k)-f|v^^(x.-e2h,t.)  ;  o<©^,92<l 

Similarly  we  obtain  from  (5»1)  t»)>  (5«1+)  b)  and  (6,1) 

(6,15) 

2 
I^(E|,Fg)=^^j=  -  X,(x^,t.-93_k).  ^u^^3^i^e2^,t.),o<©^,©2<i; 

where 

(6.16) 

=  i^  (u^,  vj)  f-3   ^   QJyj  ^  (i.e^)^(    t.) 

b.  .  =  (L^  ("J,  vJ)  [  -j  ^  ^j^j  ^  (l-0j)v(x.,  t.) 
On  the  boundary  of  the  mesh  {$,'})   we  have 


(6.17) 


E?  =  P?  =  0   .    o  <  i  <  N 


^o  =  ^o  =  ^N  =  °  '      o  <  j  <  J 


Thus  the  error  satisfies  a  system  of  equations  formally 
identical  to  (6.3),   If  condition  (6,2)  is  satisfied  then 
the  bounds  (6.12)  and  (6.13)  apply  to  |F^|  and  |E^|  re- 
spectively, where  now 


•  ■  •     I 


(6.18) 


A  =  max 


B  =  max 
ij 


4=k/2max(lv^^(x^,tj.)-^^v^x(x^,t^)|,|u^^(x^,tj)+k-^u^^i,tj)|) 

Here  the  arguments  are  understood  to  be  those  of  (6.1i4-),  (6.15) 
and  (6.16).   For  any  value  of  t,  say  t  =  T,  we  may  consider  a 
sequence  of  meshes  (5»3)  such  that  k  — *■  o,  while  Jk  =  T  and 
ti.  <  1  in  each  mesh.   Then  the  bounds  (6.12)  and  (6.I3)  hold 
in  each  mesh  but,  from  (6. 18),  4  — >o.   Thus  the  error 
approaches  zero  and  hence  the  niimerical  solution  uniformly 
converges  to  the  solution  of  the  partial  differential  equations 
(5.1)  in  R^  . 

For  the  difference  scheme  mentioned  at  the  end  of 
Section  5  we  must  add  to  (6.2)  the  condition 


(6.19) 


kA  <  1 


Then  by  the  above  argioments  it  can  be  shown  that  the  identical 
bovmds  (6,12)  and  (6.I3)  hold  for  the  error,  and  hence  the 
uniform  convergence  of  this  scheme  is  also  established. 
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